close; clear; clc;
Nt = 50; Nl = 500; 
nl = 100;
Fai = pi / 2;
u = 4*pi*10^(-7); p1 = 100;
i = sqrt(-1);
for it = 1:Nt
    w(it) = 10 ^ (it/10);
    k1(it) = sqrt(-i*w(it)*u/p1);
    for il = 1:Nl
        r(il) = il * nl;
        PE(il) = 1.0 * r(il);
        Ex(it,il) = PE(il) * p1 / (2*pi*r(il)^3) * (3*(cos(Fai))^2 - 2 + (1+k1(it)*r(il))*exp(-k1(it)*r(il)));
        ExA(it,il) = abs(Ex(it,il));
        ExB(it,il) = angle(Ex(it,il));
        Hz(it,il) = -PE(il) / (2*pi*r(il)^4*k1(it)^2) * sin(Fai)*((3 + 3*k1(it)*r(il) + k1(it)^2*r(il)^2) * exp(-k1(it)*r(il)) -3);
        HzA(it,il) = abs(Hz(it,il));
        HzB(it,il) = angle(Hz(it,il));
    end
    
end




figure(1)
loglog(w(:),ExA(:,10))

figure(2)
loglog(w(:),HzA(:,10))